#ifndef am_RHInterp_H
#define am_RHInterp_H



#include <string>
#include <tbb\concurrent_vector.h>

#include <math.h>

#include "..\AuroraTypes.hpp"
#include <borealis\Raster.hpp>
#include <borealis\AsciiRaster.hpp>
#include <borealis\Interpolation.hpp>

namespace Aurora
{

/*
Class: am_RHInterp
	Interpolates relative humidity over a spatial domain.
	Uses variables that are with respect to a plane surface of water.
	Note: Point value RH values of 0.0 is not allowed and breaks the underlying equations. Resultins in log(0) -> inf


In Variables:
		Raster<float>	DEM variable
		Raster<float>	Gridded temperature variable
		string			Relative Humidity header value from the met files
		string			Temperature header value from the met files
		string			Date header value from the met files
		string			Interpolation method
							- IDW => Inverse Distance Weighting
							- TPST => Thin plate spline with tension = 0.1
		string			Level in m to raise the RH to instead of that specified by the DEM. Good if you need to raise it to level X for cloud fraction calcs, etc.
							- "DEM" specifies to follow DEM
							- "xxx.xx" specifies the level in m

Out Variables:
		string			Name of gridded RH

Reference:
	Liston, Glen E., and Kelly Elder. 2006. A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). Journal of hydrometeorology 7: 217-234.

*/
class am_RHInterp : public AuroraTypes::ModuleBase
{
	//in vars

	int DEM;
	int TEMPGRID;
	int RH;
	int TEMP;
	int DATE;
	int INTERP;
	int LEVEL;

	// out vars
	int OUTVAR;
	std::string m_interpMethod;
public:
	am_RHInterp(std::string n)
	{
		 DEM		=0;
		 TEMPGRID	=1;
		 RH			=2;
		 TEMP		=3;
		 DATE		=4;
		 INTERP		=5;
		 LEVEL		=6;

		 OUTVAR		=0;

		 name = n;
	}


	class RaiseRHVisitor: public Borealis::RasterVisitor<float>
	{
	private:
			boost::shared_ptr<Borealis::Raster<float> > m_dem;
			boost::shared_ptr<Borealis::Raster<float> >	m_temp;
			int m_month;
			float m_level;
	public:
		RaiseRHVisitor(int month, float level, boost::shared_ptr<Borealis::Raster<float> > dem, boost::shared_ptr<Borealis::Raster<float> > tempGrid)
			: m_dem(dem), m_month(month), m_temp(tempGrid), m_level(level) {}


		void operator()(float* cell, unsigned int row, unsigned int col)
		{
			//use water for the moment as per Liston, Elder
			float a = 611.21F;
			float b = 17.502F;
			float c = 240.97F;
			float level = m_level;

		
			float lambda = 0.0F;

			//dew point lapse rate
			switch(m_month) //these are per 1m
			{
			case 1:
				lambda=0.00041F;
				break;
			case 2:
				lambda=0.00042F;
				break;
			case 3:
				lambda=0.00040F;
				break;
			case 4:
				lambda=0.00039F;
				break;
			case 5:
				lambda=0.00038F;
				break;
			case 6:
				lambda=0.00036F;
				break;
			case 7:
				lambda=0.00033F;
				break;
			case 8:
				lambda=0.00033F;
				break;
			case 9:
				lambda=0.00036F;
				break;
			case 10:
				lambda=0.00037F;
				break;
			case 11:
				lambda=0.00040F;
				break;
			case 12:
				lambda=0.00040F;
				break;
			}

			float dewPointLapseRate = lambda * c/b;

			float Td = 0.0F;
			
			if(level != -1.0F)
				Td =*cell-(-dewPointLapseRate)*(0-m_dem->GetRasterAt(row,col)); //going up
			else
				Td =*cell-(-dewPointLapseRate)*(0-level); //going up
			
			float temp = m_temp->GetRasterAt(row,col);

			float e = a * exp( (b*Td) / (c+Td));
			float es = a * exp( (b*temp) / (c+temp)); 

			//RH value replaces Tdew value
			*cell = 100*e/es;

			//not RH > 100%
			if(*cell > 100)
				*cell = 100;


		}
	};

	void run(AuroraTypes::StationList* stations,AuroraTypes::VariableList* variables, AuroraTypes::InVariables& inVars, AuroraTypes::OutVariables& outVars)
	{
		boost::shared_ptr<Borealis::Raster<float> > m_dem;
		boost::shared_ptr<Borealis::AsciiRaster> m_RH ;
		boost::shared_ptr<Borealis::Raster<float> > m_temp ;
		//bail if we don't have any met stations
		if(stations->size() == 0)
			throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) +"No met data");

		try
		{
			m_dem = GetRaster(variables,inVars[DEM]);
			m_temp = GetRaster(variables, inVars[TEMPGRID]);
		}
		catch (std::runtime_error e)
		{
			throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) +" - " + std::string(e.what()));

		}


		m_RH.reset(new Borealis::AsciiRaster(m_dem->GetRow(),m_dem->GetCol()));
		m_RH->SetXllCorner(m_dem->GetXllCorner());
		m_RH->SetYllCorner(m_dem->GetYllCorner());
		m_RH->SetCellSize(m_dem->GetCellSize());

		//for each station lower the temperature to a common level z
		//build a new point and save it for the interpolation scheme
		tbb::concurrent_vector<Borealis::RasterPoint<float> > points;

		m_interpMethod = inVars[INTERP];

		//get the date from the first station
		//we assume at this point that all the stations are the same length, have the same time step, etc
		int month = (stations->at(0))->Now().Month();

		//correct the reference temperatures
		for(unsigned int i = 0; i < stations->size();i++)
		{
			float rh = stations->at(i)->Now().Get<float>(inVars[RH]);
			float temp = stations->at(i)->Now().Get<float>(inVars[TEMP]);
			unsigned int col = stations->at(i)->GetX();
			unsigned int row = stations->at(i)->GetY();

			float newTd=0.0F;

			

			//use water for the moment as per Liston, Elder
			float a = 611.21F;
			float b = 17.502F;
			float c = 240.97F;
			
			//solve RH ~= 100* e/es for e
			float es = a*exp( (b*temp)/(c+temp) );
			float e = rh / 100  * es;

			//we must not allow rh to be == 0.0, other wise the above eqn fails
			if(rh < 0.0)
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + " bad RH value. RH = 0.0");	

			//calculate dewpoint temperature
			float dewPointTemp = (c*log(e/a)) / (b-log(e/a));

			float lambda = 0.0F;

			//dew point lapse rate
			switch(month) //these are per 1000m
			{
			case 1:
				lambda=0.00041F;
				break;
			case 2:
				lambda=0.00042F;
				break;
			case 3:
				lambda=0.00040F;
				break;
			case 4:
				lambda=0.00039F;
				break;
			case 5:
				lambda=0.00038F;
				break;
			case 6:
				lambda=0.00036F;
				break;
			case 7:
				lambda=0.00033F;
				break;
			case 8:
				lambda=0.00033F;
				break;
			case 9:
				lambda=0.00036F;
				break;
			case 10:
				lambda=0.00037F;
				break;
			case 11:
				lambda=0.00040F;
				break;
			case 12:
				lambda=0.00040F;
				break;
			default:
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + "unknown month for dewpoint lapse rate calculation");				
			}

			float dewPointLapseRate = lambda * c/b;

			
			newTd =  dewPointTemp - dewPointLapseRate*(0.0 - stations->at(i)->GetElevation());

			Borealis::RasterPoint<float> point(row,col,newTd);
			points.push_back( point );
		}

		std::string l = inVars[LEVEL];
		float level = -1.0F;

		if(l != "DEM")
		{
			try
			{
				level = boost::lexical_cast<float>(l);
			}
			catch ( ... )
			{
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + " bad RH level");				
			}
			
		}


		if(m_interpMethod == "IDW")
		{
			RaiseRHVisitor* visitor = new RaiseRHVisitor(month, level, m_dem, m_temp );
			Borealis::Interpolation::InverseDistanceWeighting(points,  m_RH, visitor);
			delete visitor;
		}
		else if(m_interpMethod == "TPST")
		{
			RaiseRHVisitor* visitor = new RaiseRHVisitor(month, level, m_dem, m_temp );
			Borealis::Interpolation::ThinPlateSplineT(points,  m_RH, visitor);
			delete visitor;
		}
		else
		{
			throw std::runtime_error("Unknown interpolation method");
		}

		AuroraTypes::VariableList::accessor a2;
		variables->insert(a2,outVars[OUTVAR]);		
		a2->second = m_RH;

	}


};
}

#endif